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Abstract 

We study the transverse cracking of the 90° ply in [0/90] s cross-ply laminates 
by means of a discrete element method. To model the 90° ply a two-dimensional 
triangular lattice of springs is constructed where nodes of the lattice model fibers, 
and springs with random breaking thresholds represent the disordered matrix ma- 
terial in between. The spring-lattice is coupled by interface springs to two rigid 
bars which represent the two 0° plies in the model, which could be sublaminate as 
well. Molecular dynamics simulation is used to follow the time evolution of the 
model system. It was found that under gradual loading of the specimen, after some 
distributed cracking, segmentation cracks occur in the 90° ply which then develop 
into a saturated state where the ply cannot support additional load. The stress 
distribution between two neighboring segmentation cracks was determined, fur- 
thermore, the dependence of the microstructure of damage on the ply's thickness 
was also studied. To give a quantitative characterization of stiffness degradation, 
the Young modulus of the system is monitored as a function of the density of seg- 
mentation cracks. The results of the simulations are in satisfactory agreement with 
experimental findings and with the results of analytic calculations. 
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1 Introduction 

Transverse cracks of 90° plies in [0/90]s cross-ply laminates [?] appear early under 
tensile loading and continue forming until reaching a saturation state. The understand- 
ing of multiple transverse cracking is not only a phenomenon of scientific interest but 
also of practical importance. Transverse cracks can induce other damage mechanisms 
such as micro-delaminations followed by micro-buckling [?]. Failure generally occurs 
by the growth of nucleated or intrinsic defects on length scales much smaller than the 
sample size or representative volume element size used for the description of trans- 
verse cracking. One goal in materials design is the production of materials with flaw- 
insensitive or damage tolerant behavior, in a way that small flaws propagate in a stable 
manner up to global material failure [?]. This necessitates the activation of interacting 
stress releasing mechanisms like crack deflections out of the preferred crack growth 
plane, micro-delaminations, induced microcracks or crack bridging. In any case the 
crack no longer propagates inside a homogeneous continuum, and inhomogeneity or 
disorder on certain length scales becomes a necessary element for the failure predic- 
tion of the composite. Models for the fracture simulation should allow the study of the 
growth of small cracks up to macroscopic size. 

For the last 20 years several analytical solutions were proposed for the calculation 
of the density of micro cracks under mechanical and thermal loading, using a variety 
of different approaches for stress analysis and different failure models [?, ?, ?]. These 
models were combined with interacting damage mechanisms within a micro-mechanics 
of damage analysis to take micro-delaminations into account [?]. For woven fabric sim- 
ple models were derived from cross-ply laminates, showing the same types of damage 
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also found in cross-ply laminates. Even though analytic models are in good agreement 
with experimental findings, the description of the evolution of transverse cracks from 
the early nucleation of micro-cracks over coalescence and crack growth lies beyond the 
scope of an analytic approach. 

Numerical simulations for the failure of the transverse ply have been performed at 
different length scales ranging from the microscopic scale, modeling single fibers [?, ?] 
and models for small fiber clusters [?, ?, ?] up to the meso scale studying homogenized 
plies [?] with a variety of model approaches. Either symmetric fiber arrays built of unit 
cells consisting of a part of a fiber and matrix material are calculated using finite ele- 
ment methods [?, ?] or small clusters of fibers with microstructural disorder, embedded 
in matrix material are modeled directly, employing fictitious or discrete crack models 
[?]. Unfortunately, these simulations are very time consuming and the results have 
limited meaning for the simulation of multiple cracking in thicker transverse plies. In 
addition, they do not allow predictions about size scaling or the possibility to model 
the dynamic fracture process itself. 

In the present paper we study the transverse cracking of the 90° ply in a [0/90] s 
cross-ply laminate by means of a discrete element method. In order to model the 90° 
ply, a two-dimensional triangular lattice of springs is constructed. The nodes of the lat- 
tice model fibers perpendicular to the plane of the lattice, and springs represent the ma- 
trix material in between. The spring-lattice is coupled by interface springs to two rigid 
bars which represent the two 0° plies in the model, but could be adjacent sublaminates 
as well. Disorder is introduced by assigning randomly distributed breaking thresholds 
to the springs, i.e. a spring breaks if the load on it exceeds its breaking threshold. To 
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simulate transverse cracking under uniaxial loading of the composite, external load is 
imposed on the triangular lattice by horizontally displacing the two edges and the up- 
per and lower interface elements of the lattice. This corresponds to a load imposed by 
coupled flexural rigid bars along the load direction. The time evolution of the system is 
followed by solving the equation of motion of the nodes (molecular dynamics). With 
the discrete element method used in this study, relatively large system sizes can be 
handled, enabling simulation of the fracture processes with crack-crack interactions of 
thousands of cracks. Multiple cracking has been studied before, using random spring 
networks with springs of random distribution of strain failure [?, ?, ?, ?], but rarely in 
the field of fiber composites research. This approach is particularly suited for studies 
on dynamic instability in crack propagation, the collective behavior of many interact- 
ing cracks and size effects of multiple transverse cracking. Therefore, the discrete 
element method is practical for studying general features of the statics and dynamics 
of fracturing like the crack morphology, global fracture patterns due to the collective 
behavior of many interacting cracks as well as the dynamic instability in the propa- 
gation of cracks. Our investigation is focused on the process of damaging, stiffness 
degradation, micostructure of damage, and furthermore, on the relative importance of 
damage mechanisms like segmentation and delamination in the degradation process. 

The paper is organized as follows: Section 2 gives a detailed description of the 
model construction. Tests of the model and the numerical results on transverse crack- 
ing obtained by simulations are presented in Section 3. An analytic approach worked 
out in the literature to the transverse cracking problem is summarized in Section 4. 
The simulational results are confronted with the analytic calculations and experimental 
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results at the presentation of the analytic methods and also in the discussion part. 

2 The discrete element model 

In order to study the formation of segmentation cracks, a two dimensional triangular 
spring network model of composites containing undeformable cells, is worked out. 
Molecular Dynamics (MD) simulation is used to follow the motion of each cell by 
solving Newton's equations of motion. In the present study we use a 4 t/l order Gear 
Predictor Corrector scheme [?]. A general overview of MD simulations applied to 
composite materials can be found in [?]. 

With this method already small cracks are sharply defined, with the possibility 
to simulate simultaneously a conglomerate of cracks within rather large lattice sizes. 
The fundamental advantage of the lattice model used in this investigation is due to 
its simplicity, giving direct access and possible physical interpretation to each step of 
the algorithm. Consequently, one can modify the rules for features of interest, like 
the characteristic properties of size, strength or force in a rather straightforward and 
transparent way. 

Despite of the advantages of lattice models, one has to be aware of the limitations of 
discrete network models for the simulation of failure. One type of questions deals with 
the abstraction of the underlying material regions with discrete elements concerning the 
attributes of the elements and their organization in a network as well as the boundary 
conditions and numerical approaches utilized. But these questions are independent 
from limitations imposed by the nature of fracture simulations with discrete elements. 
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Discreteness implies, that the stresses in an element correspond to an average of the 
local stresses over the volume represented by this element. It is, therefore, impossible 
to model cases with stress transfer from failed material regions to the neighboring ones 
on a length scale smaller than the lattice spacing. 

The model is composed in three major steps, namely, a) the implementation of the 
microscopic structure of the solid, b) the determination of the constitutive behavior, 
and finally c) the breaking of the solid. 

(a) Microstructure The model is composed of circular cells of identical radius 77 
with their centers located on the nodes of a regular triangular lattice built out of equi- 
lateral triangles with s being the characteristic lattice spacing as illustrated in Fig. 1. 
The cells represent the cross sections of fibers. The fiber radius 77 depends on the fiber 
volume fraction Vf in the transverse ply, and is adjusted following the expression 



The cell centers are connected by linear elastic springs with stiffness E s , represent- 
ing the material region with fiber, fiber matrix interface and embedding matrix mate- 
rial. For simplification purposes, the springs behave in a time independent way, but 
time dependent spring properties are easy to introduce. In our simulation, the cells are 
the smallest particles interacting elastically with each other. We use a two-dimensional 
simulation with cells of unit thickness, allowing only motion in the observation plane, 
with the two degrees of freedom being the two coordinates of the position of the cell 
center of mass. The utilization of a regular lattice is a clear neglect of the topolog- 



r f = 




< v f < v fmax = 0.906. 



(1) 
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ical disorder, symptomatic for realistic fiber composite systems. Since we have no 
particular basis for a topological disorder, we make a computational simplification by 
assuming the disorder in spring properties to adequately account for all the relevant dis- 
order present in the material. Disorder is given to the model by a Weibull distribution 
of breaking thresholds Fd at the beginning of the simulation in the form 



and is kept constant during the fracture process (quenched disorder). The Weibull 
modulus m controls the degree of disorder in the distribution and it is usually chosen 
in the range 1.5 < m < 20, experimentally found to describe a variety of materials. F 
represents the characteristic strength of the material. The top and bottom row of cells, 
as well as the left and right rows are connected to flexural rigid beams via interface 
springs, representing the 0° — 90° interface of the composite. In the following nx 
denotes the number of cells along the loading direction, and ny along the perpendicular 
one. 

(b ) Constitutive behavior The fibers are assumed to be rigid undeformable cylindrical 
bodies with the possibility to overlap when they have contact. This overlap represents 
to some extent the local deformation of the fibers (soft particle approximation [?]). A 
contact law is used to estimate the contact force F~ between two circular cells i and j 
as a function of the overlap area A (see Fig. 16). The overlap area A can be expressed 




(2) 
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in terms of the distance r, j of the centers as 




(3) 



and the repulsive contact force between two cells is given by 




EftA 




(4) 



where L c is a characteristic length chosen to be equal to the lattice spacing s, Eft 
denotes the elastic modulus of a fiber perpendicular to the fiber direction, and riy is 
the unit vector pointing from particle i to j. The contact law inside the system is 
only applied between two circular particles after the breaking of the connecting spring, 
preventing penetration of broken parts into each other. It is a drawback of this method 
that the overlap force (^) does not have a Hamiltonian formulation which implies that 
no energy can be associated to the deformation modeled by the overlap. The force Fi on 
the cell i is the residual force of all attractive and repulsive forces on the cell calculated 
from their interaction with neighboring cells either by contact or by the springs. In 
principle, stress can be applied to the system either by changing the elastic modulus, 
the initial length of the springs or by applying a force on the interface elements. In 
the present paper for the purpose of the study of transverse ply cracking, mechanical 
load is introduced on the transverse ply by straining the rigid upper and lower beams 
in x-direction and moving the left and right beam, respectively (see also Fig. la). 
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(c) Breaking of the solid Stresses inside the system can be released either by the 
formation of surfaces (cracks) inside the transverse ply or by fracture of the interface 
elements on the boundary. In the framework of Discrete Element Methods, the com- 
plicated crack-crack interaction is naturally taken into account. If the total force on a 
spring F exceeds its breaking threshold F4 (damage threshold), its stiffness is abruptly 
reduced to zero, resulting in a load redistribution to neighboring springs in the next 
iteration steps. If after some iterations the neighboring springs exceed their threshold 
value due to the additional load, they fail too, giving rise to crack growth. Like in the 
bulk material, the interface springs break once their maximum load is exceeded. Their 
failure is interpreted as micro delamination. The breaking characteristics of the inter- 
face elements are in principle independent of those of the bulk springs. In the model 
this is captured by assigning damage thresholds to the interface springs using a Weibull 
distribution but with parameters different from the ones used for the bulk. 

The time evolution of the particle assembly is followed by solving the equation of 
motion of the nodes {i.e. transverse fibers) 

mf'i =Fi, i = 1, ... ,N (5) 

where N denotes the number of nodes (fibers of the transverse ply), and Fj is the 
total force acting on node i. A 4 th order Predictor Corrector algorithm is used in the 
simulations to solve numerically the second order differential equation system Eq. (|j). 
After each integration step the breaking condition is evaluated for all the intact springs. 
Those springs which fulfill this condition are removed from the further calculations 
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(spring breaking). 



3 Computer simulations 

Before applying the model to study transverse cracking in fiber composites a variety 
of simulations have been performed in order to test the behavior of the model and to 
make parameter identification with respect to experiments. Therefore, reasonable char- 
acteristic values for the lattice spacing L c , spring stiffness E s and strength <r c , contact 
modulus Ef t and time have to be chosen. The Young modulus and Poisson ratio of the 
model system was measured numerically by simulating uniaxial loading of a rectangu- 
lar sample when springs were not allowed to break. To avoid the disturbing effect of 
elastic waves that can be induced by the loading in our dynamical model, the numeri- 
cal experiment is performed in the way developed in [?]: The two opposite boundaries 
of the solid are moved with zero initial velocity but non-zero acceleration. When a 
certain velocity is reached the acceleration is set to zero and the velocity of the bound- 
aries is kept fixed, which ensures a constant strain rate loading. With the help of this 
slow loading the vibrations of the solid can be reduced drastically compared to the case 
when the boundaries start to move with non-zero initial velocity. An additional way to 
suppress artificial vibrations in dynamical models, used in this study, is to introduce a 
small dissipation (friction or damping) for the springs. This dissipation has to be small 
enough not to affect the quasi-static results. In order to measure the macroscopic elastic 
behavior of the model solid, the total force needed to maintain the external constraint 
is monitored from which the stress a can be determined as a function of strain e. A 
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representative example of the macroscopic elastic behavior is shown in Fig. 2a, where 
a good agreement of the numerical and analytical behavior can be observed. 

As a further test of the model, at a certain load level two parallel cracks were intro- 
duced in the lattice perpendicular to the loading direction spanning through the whole 
cross section of the transverse ply. Keeping fixed the external load, the distribution of 
the stress components along the interface between the transverse and the longitudinal 
plies was calculated. The distribution of a xx is presented in Fig. 26. 

3.1 Crack formation under gradual loading 

The discrete element method gives the possibility to monitor the development of the 
microscopic damage in the specimen. Simulations showed that in accordance with ex- 
periments the damage process of the transverse ply is composed of two distinct parts. 
At the beginning of loading the weakest springs break in an uncorrected fashion when 
they get over-stressed, generating microcracks in the specimen. This primary uncor- 
rected microcrack nucleation is dominated by the disorder of the system introduced 
for the damage thresholds of the springs. The microcracks substantially change the 
local stress distribution in the ply leading to stress concentrations around failed re- 
gions, which gives rise to correlated growth of cracks under further loading, see Fig. 
3a. Crack growth occurs perpendicular to the loading direction. Growing cracks or 
cracks formed by the coalescence of growing cracks can span the entire thickness of 
the transverse ply, resulting in segmentation, i.e. break-up of the ply into pieces, as is 
illustrated in Fig. 36. Reaching the ply interface the crack stops without the possibility 
of penetration into the longitudinal ply in the model, as was observed in most of the 
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experiments, but with the possibility of deflection at the ply interface. The simulated 
fracture process is similar to observations described in the literature [?]. 

Further segmentation cracks mainly form between existing cracks until the crack 
density is saturated due to occurring micro delaminations. Micro delaminations usually 
start to occur when the crack density has already reached high values, depending on the 
thickness of the transverse ply. Fig. 3d,e show that as a result of the formation of cracks 
spanning the entire thickness, the ply gets segmented into pieces, with delamination 
zones extending along the interfaces. 

In order to obtain a quantitative characterization of the microstructure of damage, 
we analysed the spatial distribution of microcracks, i.e. the distribution of broken 
springs. A crack is identified as a connected set of broken springs in the triangular 
lattice taking into account solely nearest neighbor connections. A crack is considered 
to be a segmentation crack if it spans from one side of the transverse ply to the other 
one. Since cracks forming along the interface of the plies complicate the identification 
of segmentation cracks in the framework of our algorithm, we use now the inner part of 
the specimen for the cluster identification. The method is illustrated in Fig. 4 for two 
different snapshots of the loading process. Fig. 4a and Fig. 4c contain all the broken 
springs which are present at the given time, while in Fig. 46 and in Fig. Ad only the seg- 
mentation cracks are presented which are identified by our analysis. It can be observed 
that the number of microcracks and also the number of segmentation cracks increases 
with time and the segmentation cracks have a more or less quasi-periodic spacing. The 
formation of the segmentation cracks is very fast and often shows stable crack growth 
only during the crack initiation. After the formation of the first segmentation crack 
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the specimen is strained further, keeping the velocity of the boundary cells constant, 
till enough energy is built up to activate the formation of the next segmentation cracks 
approximately midway between existing segmentation cracks. 

The total number of microcracks divided by the total number of bonds in the system 
N cracks /Nb on d s and the number of segmentation cracks N seg are presented in Fig. 5 
as a function of strain E s e/<r c for several different values of the Weibull-modulus m. 

It can be observed that cracking initiates at a finite strain value e in , called damage 
initiation strain, and starting from this point N cr acks monotonically increases during 
the entire loading process. Segmentation first occurs at a strain larger than £,„ after the 
number of microcracks reached a certain value. The number of segmentation cracks 
N seg also increases, however, it gets saturated for larger strains in accordance with 
experimental observations. 

Fig. 5 provides also an insight into the role of disorder in the damage process. 
Increasing the value of m, i.e. making the Weibull distribution narrower, which implies 
that the system gets more and more brittle, the damage initiation strain e in and the 
corresponding strain value at which the segmentation starts increase, furthermore, the 
saturation value of the segmentation cracks also increases for more brittle systems. The 
slight decrease of N seg for wider breaking threshold distributions is an artifact of the 
cluster analysis that identifies segmentation cracks connected via large delamination 
zones as one cluster and therefore as one segmentation. 
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3.2 Transverse ply stiffness degradation 

All the observed forms of damage lead to a reduction of the total load portion, the 
transverse ply shares with the longitudinal plies. The ply stiffness is a macroscopic 
property, but the numerical model is of micro structural nature. Therefore, the different 
effects of damage on the global stiffness have to be obtained by calculating an effective 
stiffness of the transverse ply. In this study we calculate the effective Young's modulus 
E e f f by cumulating the energy density e pot stored in the spring elements as 

Eeff = ■ (6) 

This assumption is justified due to the way of loading, resulting in mainly elongated 
springs, and it is also ensured that only a small portion of the cell-cell contacts are clos- 
ing in the late stage of the simulations. Fig. 6 shows stiffness degradations with varying 
Weibull modulus and with cracking artificially switched off after a certain number of 
iteration steps in the simulations. It can be observed that the more brittle the system is 
(i.e. the higher the Weibull modulus m is), the longer is the horizontal plateau that pre- 
cedes the decreasing part of the effective young modulus of ply group 1 E^j, which 
starts at the crack initiation strain. For high strains all the curves with different m values 
tend to a finite but very small stiffness which is reached when the cracking stops due to 
saturation. When cracking is artificially switched off after a certain time, the effective 
stiffness remains constant due to the constant amount of damage as expected. This 
corresponds directly to a change in the secant stiffness, useful for continuum damage 
models. 
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3.3 Micro-cracking in laminates with varying cross-ply thickness 

The damage evolution in the transverse ply of [0, 90„] s laminates depends not only on 
the disorder, but also on the thickness of the transverse ply. It is observed experimen- 
tally that for thicker transverse plies, the saturation crack density is smaller than for 
thin plies, with delamination starting earlier during the loading. Both effects can be 
observed qualitatively in Fig. 7, two compare snapshots of the damage process with 
different thicknesses see also Fig. 3a — e. 

4 Analytical model 

Popular analytical models for the composite damage analysis are continuum damage 
mechanical (CDM) or micro-mechanics of damage (MMD) models. In CDM, the me- 
chanical properties are expressed as functions of damage, with the aim of finding the 
constitutive relations between different states of stress, strain and damage. The rep- 
resentation of all kinds of damage within one damage tensor, expressing the state of 
damage that may include a combination of matrix cracking, fiber breakage, micro de- 
lamination and ply delamination, shows the dependency and limitation of CDM from 
experimental and/or theoretical input. The main disadvantage of a CDM approach for 
this study is believed to be, that CDM does not use crack propagation theories based 
on fracture mechanics, which is necessary for a damage propagation description. 

MMD analysis is based on the individual prediction of the initiation and propaga- 
tion of the various types of damage, including possible interactions between all identi- 
fied relevant damage mechanisms. Such analysis consist of two steps, a stress analysis, 
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and a damage analysis. The stress analysis has to be undertaken in the presence of ob- 
served damage. The stress analysis is completely independent from the postulation of 
a failure criterion for the initiation and evolution of damage, but of course not from the 
damage itself. An overview of MMD analysis for matrix micro cracking can be found 
in [?, ?]. Since the analytic solution is in good agreement with experimental results, 
we use it for comparison with results of this study [?, ?]. 

Stress analysis The stress field in presence of micro-cracking damage is calculated 
for a unit cell with cracks located at each end (Fig. 8). The entire specimen is built 
by a sequence of unit cells. A two-dimensional analysis for the x — z— plane derived 
by Hashin [?, ?], using a variational mechanics analysis is employed. This way, the 
characteristic violation of the boundary condition caused by a non-zero interfacial shear 
stress at the surface of the ply crack, which is symptomatic for an one-dimensional 
analysis can be avoided. Following [?] we assume that the x-axis tensile stress in each 
ply group is independent of z, that ply cracks span the entire ply thickness and that 
delaminations are always symmetric [?]. 

The x-axis tensile stress of a micro-cracked laminate can be derived with force 
balance 

<&> = EW/E%* - 1>{x), 4 2 i = 4 2) /£ c % + (7) 

with the undetermined function ij){x) and A = t 2 /ti expressing stress perturbations 
caused by the segmentation cracks. is again the x-direction young's modulus of 
ply group i and denotes global laminate properties. The integration of the stress equi- 
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librium equations with the unit cell boundary conditions gives the shear and transverse 
stresses in the unit cell in terms of ij)(x). The function ip(x) that minimizes the comple- 
mentary energy gives the best approximation to the micro-cracked cross-ply laminate 
stress state. The solution of the Euler equation for finding tp(x) can be found in detail 
in [?]. Fig. 9 shows the stress distributions inside a unit cell for interfacial values of 
a xz and <j zz as well as a xx for longitudinal and transverse ply, with the comparison 
of crix to the stress calculation from our numeric simulation (compare Fig. 9a). A 
reasonable agreement of the analytic solution and of the corresponding simulation can 
be observed. 

The composite properties for the laminate have been calculated by means of micro- 
mechanics from the bulk properties for the fiber and matrix constituents [?, ?]. For 
this work we renounce a detailed tabulation of bulk and composite material properties, 
since emphasis is given on a dimensionless presentation of results. Nevertheless, the 
characteristic length of the numerical model is scaled with the fiber volume fraction and 
an average fiber diameter of 7/j.m, whenever it was necessary for comparative reasons. 

Failure model The failure criterion used in the analytical model, is an energy cri- 
terion [?], allowing a formation of a new crack, whenever the strain energy release 
rate (ERR) G m associated to the formation of a new micro crack exceeds some criti- 
cal value G mc , also called intralaminar fracture toughness or micro-cracking fracture 
toughness. The ERR is also calculated with Hashin's two-dimensional, variational me- 
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chanics analysis [?]. 



G m = <4o f l V C l (^4 - ^2 + 2VC1C3. (8) 



Cl - T7TT + TT72T' C2 ~ T7TT A + o ~ 



(9) 



1 A 3 1 A 
— (15A 2 + 20A + 8) + m , Ci = ^- + ^ 



E 1 ^ . etc, are the material properties of the plies in conventional notation. 



Delamination Microcrack-induced delamination is a damage mechanism, that com- 
petes with the formation of new segmentation cracks in the transverse ply. Experi- 
mental and numerical studies show, that ply crack-induced delamination does not start 
before a critical crack density is reached. 

The two-dimensional, variational mechanics analysis of [?] described above can be 
extended to account for delaminations emanating from ply crack tips. The assumption 
that trix is independent of z, indicates identical delamination lengths emanating from 
the top and bottom crack tips. This leads to a simple stress state for the delaminated 
regions: 

4^=0, 42 = <&> = (), (id 

4 2 i=4 2 o } +^f = ^o, 4 2 ;=4 2 ;=o. (12) 
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In the region between the delaminations, the stress distribution is equal to that existing 
between two microcracks with the dimensionless half spacing p — 5 with S = (dl + 
d2)/2t\ and p = a/t\. The ERR for the initiation of delaminations in the unit cell is 
given by: 



x(p) is a function that corresponds to the excess strain energy in a unit cell of damage 
caused by the presence of ply cracks [?, ?]. Assuming that the critical energy release 
rate for delamination G mc — Gdc, we compare the dimensionless ERR as a function 
of the crack density to find the crack density for the initiation of delaminations [?]. 
In the analytic framework, with the onset of delamination, no further segmentation 
takes place, so the crack density for delamination onset is equal to the crack density of 
the characteristic damage state (CDS). In Fig. 9b the analytic values for the CDS as a 
function of transverse ply thickness are compared to our simulation results and to some 
limited experimental data. Since our numerical model is more realistic in the sense that 
it captures more details of the microscopic damage mechanisms, a satisfactory agree- 
ment can be stated between the simulation and the experimental findings, however, the 
analytic curve falls a bit further from them. 

The reduction of the transverse ply stiffness is compared to simulation results in 
Fig. 10. It is obvious, that the analytic model without delamination has no segmen- 
tation crack saturation, leading to a higher crack density. Numerical results for two 
different transverse ply thicknesses show qualitative agreement, only with the stiffness 




(13) 
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reduction being zero at the saturation crack density. 

5 Discussion 

We introduced a disordered spring network model to study the transverse cracking of 
the 90° ply in [0/90 n ]s cross-ply laminates. The main advantage of our modeling is 
that it naturally accounts for the complicated local stress fields formed around failed re- 
gions in materials, furthermore, it captures the gradual activation of the relevant failure 
mechanisms and their interactions during the fracture process. 

We have demonstrated that our discrete element model provides a deep insight into 
the damage process occurring under gradual loading of cross ply laminates. Quanti- 
tative results have been obtained on the microstructure of damage, on the local stress 
distribution between cracks and on the degradation of the ply's stiffness. 

The results obtained by numerical simulations have been confronted with experi- 
mental findings [?] and also with the analytic approach of Nairn and Hu. Good agree- 
ment was found between the numerical, analytical and experimental results. How- 
ever, the numerical simulations proved to be more realistic than the simple analytic 
approaches, due to the more detailed description of microscopic damage mechanisms 
as built in the discrete model. 

Despite of the realistic model construction our study makes some technical simpli- 
fications which have to be overcome in the future. For instance, due to the difficulties 
in the identification of the fracture properties of interface bonds between plies, no sim- 
ulations have been performed varying the relative breaking thresholds of interface and 
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bulk springs. To clarify this problem a detailed parameter study will be carried out in a 
forthcoming publication. 

The capabilities of our model are not limited to uniaxial loading, it can be applied 
to study damage processes of the transverse ply in cross-ply laminates under varying 
loading cases including also the thermal degradation process occurring during the pro- 
duction process of fiber composites. Studies in this direction are in progress. 
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Figure Legends 



Figure 1 : Microstructure of the model (a) and cells in contact (b) 

Figure 2: (a) Macroscopic constitutive behavior of the uncracked specimen in 
a dimensionless form, obtained by uniaxial loading (switching off the spring 
breaking in the simulation), (b) Stress distribution <r xx /a x0 between two cracks. 

Figure 3: Snapshots of the model system with the size nx = 800, ny = 10. (a) 
nucleated microcracks, distributed over the whole test cell, and the formation of 
the first segmentation crack. (6) Quasi periodic segmentation pattern with cor- 
responding stress distribution (c). (d) Segmentation cracks and delaminations at 
crack density saturation state, (e) The saturated state for a system of nx = 800, 
ny = 50 

Figure 4: Identification of segmentation cracks at two different stages of the 
damage process, (a) and (c) show all the microcracks which are present at a 
given state, in (6) and (d) only the identified segmentation cracks are presented. 

Figure 5: (a) The total number of microcracks N cra cks divided by the total 
number of bonds Nb on ds, and (6) the number of segmentation cracks N seg as a 
function of strain E s e/a c . Ply thickness ny was 10 cells. Smooth curves were 
obtained by averaging over six samples. 

Figure 6: Stiffness reduction with varying Weibull moduli (a) and for m = 4 
switching off the breaking after a certain number of iteration steps. 
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Figure 7: (a) The total number of microcracks N crac k s divided by the total 
number of bonds Nb on ds, and (6) the number of segmentation cracks N seg as a 
function of strain E s e/<r c for systems of changing thickness ny. 

Figure 8: Side view of an unit cell representing a cross ply laminate with seg- 
mentation cracks and delaminations. 

Figure 9: Stress distribution inside the unit cell for a = 2t\ and t\ = ti (a), and 
characteristic crack density, analytic results for G Cdel = G Cseg , numerical and 
experimental [?] results (b). 

Figure 10: Comparison of analytic predictions for the loss of stiffness E^/E^ 
with simulation results by our model. 
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Figure 2: (a) Macroscopic constitutive behavior of the uncracked specimen in a di- 
mensionless form, obtained by uniaxial loading (switching off the spring breaking in 
the simulation), (b) Stress distribution (J xx /<t x0 between two cracks. 
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Figure 3: Snapshots of the model system with the size nx — 800, ny = 10. (a) 
nucleated microcracks, distributed over the whole test cell, and the formation of the 
first segmentation crack, (b) Quasi periodic segmentation pattern with corresponding 
stress distribution (c). (d) Segmentation cracks and delaminations at crack density 
saturation state, (e) The saturated state for a system of nx — 800, ny = 50 
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Figure 4: Identification of segmentation cracks at two different stages of the damage 
process, (a) and (c) show all the microcracks which are present at a given state, in (b) 
and (d) only the identified segmentation cracks are presented. 
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Figure 5: (a) The total number of microcracks N crac k s divided by the total number of 
bonds Nbonds, and (b) the number of segmentation cracks N seg as a function of strain 
E s e/<r c . Ply thickness ny was 10 cells. Smooth curves were obtained by averaging 
over six samples. 



30 



F. K. Wittel, F. Kun, B. H. Kroplin, H. J. Herrmann: A 'Study of Transverse Ply 

Cracking Using a Discrete Element Method 




Figure 6: Stiffness reduction with varying Weibull moduli (a) and for m = 4 switching 
off the breaking after a certain number of iteration steps. 
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Figure 8: Side view of an unit cell representing a cross ply laminate with segmentation 
cracks and delaminations. 
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Figure 9: Stress distribution inside the unit cell for a = 2t\ and^i = t2 (a), and charac- 
teristic crack density, analytic results for G Cdcl = G Cse , numerical and experimental 
[?] results (6). 
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Figure 10: Comparison of analytic predictions for the loss of stiffness E]./E]. with 
simulation results by our model. 
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